###############################################################################
# R-Script to replicate:
# Substantive Representation of Women: Empirical Evidence
# British Journal of Political Science
# Authors: YVES KLÄY, DAVID STADELMANN, MARCO PORTMANN AND REINER EICHENBERGER
# Date: August 30, 2024
###############################################################################

#---------------------
# A) Install and Load Packages
#---------------------
install.packages("tidyverse")
install.packages("dbplyr")
install.packages("fixest")
install.packages("ggpubr")
install.packages("vtable")
install.packages("jtools")
install.packages("scales")

library(tidyverse)
library(dbplyr)
library(fixest)
library(ggpubr)
library(vtable)
library(jtools)
library(scales)

#---------------------
# B) Data and general settings
#---------------------

setwd(" - set directory where files are located - ")

source("discrete_effects.R", echo = TRUE)

parl_votes <- readRDS("femalerepresentation_data.rds")

#---------------------
# C) Functions
#---------------------
# This function calculates the discrete effects of a difference in women's "yes" vote share 
# from 0.45 to 0.55 on the decisions of female and male legislators using the DiscreteEffect function.

discrete_triple <- function(fit, lower_bound = 0.45, upper_bound = 0.55, level = 0.95)
{
  eff_male <- DiscreteEffect(fit, Probabilities = list(Prob1 = list(vox_female_vote = lower_bound),
                                                       Prob2 = list(vox_female_vote = upper_bound)),
                             FixedValues = list(female = 0),
                             Level = level
  )
  
  eff_female <- DiscreteEffect(fit, Probabilities = list(Prob1 = list(vox_female_vote = lower_bound),
                                                         Prob2 = list(vox_female_vote = upper_bound)),
                               FixedValues = list(female = 1),
                               Level = level
  )
  
  eff_diff <- DiscreteEffect(fit, Probabilities = list(Prob1 = list(vox_female_vote = lower_bound, female = 0),
                                                       Prob2 = list(vox_female_vote = upper_bound, female = 0),
                                                       
                                                       Prob3 = list(vox_female_vote = lower_bound, female = 1),
                                                       Prob4 = list(vox_female_vote = upper_bound, female = 1)),
                             LinkingFormula= ~ (Prob4 - Prob3) - (Prob2 - Prob1),
                             Level = level)
  
  rbind(eff_male, eff_female, eff_diff) %>%
    add_column(effect = c("male legislator", "female legislator", "difference"), .before = 1) %>% 
    return()
}


conf_multiplier <- qnorm(0.975)


#---------------------
# D) Descriptive Statistics
#---------------------
# Calculations for numbers used in text

# Number of MPs
n_distinct(parl_votes$mp_id)

# Number of MPs per house
parl_votes %>% 
  distinct(mp_id, council_id) %>%
  mutate(council_id = if_else(council_id==1, "UH", "LH"),
         present = 1) %>%
  pivot_wider(names_from = council_id, values_from = present, values_fill = 0) %>%
  summarize(only_UH = sum(UH==1 & LH==0),
            only_LH = sum(UH==0 & LH==1),
            both_houses = sum(UH==1 & LH==1))

# Number of referenda
n_distinct(parl_votes$ref_id)

# Number of social policy referenda
n_distinct(parl_votes$ref_id[parl_votes$top_10_social == 1])

# Share of female MPs
n_distinct(parl_votes$mp_id[parl_votes$female == 1])/n_distinct(parl_votes$mp_id)

# Share of female MPs in 1996
n_distinct(parl_votes$mp_id[parl_votes$female == 1 & parl_votes$year == 1996])/
  n_distinct(parl_votes$mp_id[parl_votes$year == 1996])

# Share of female MPs in 2022
n_distinct(parl_votes$mp_id[parl_votes$female == 1 & parl_votes$year == 2022])/
  n_distinct(parl_votes$mp_id[parl_votes$year == 2022])


# Number of referenda in which majority of women and men did not vote in the same way

# All policies
# absolute
n_distinct(parl_votes$ref_id[parl_votes$vox_gender_majority_diff == 1])
# relative
n_distinct(parl_votes$ref_id[parl_votes$vox_gender_majority_diff == 1])/
  n_distinct(parl_votes$ref_id)

# Social Policies
# absolute
n_distinct(parl_votes$ref_id[parl_votes$vox_gender_majority_diff == 1 & parl_votes$top_10_social == 1])
# relative
n_distinct(parl_votes$ref_id[parl_votes$vox_gender_majority_diff == 1 & parl_votes$top_10_social == 1])/
  n_distinct(parl_votes$ref_id[parl_votes$top_10_social == 1])

# Median of cantonal mean absolute deviation in constituency preferences across all referenda
parl_votes %>%
  distinct(ref_id, canton_id, .keep_all = TRUE) %>%
  select(yes_share_canton, ref_id, canton_id) %>%
  group_by(ref_id) %>%
  summarize(mad_ref = mad(yes_share_canton)) %>%
  summarize(med_mad_ref = median(mad_ref))

#------------
# i) Figure 1
#------------
# Differences between female and male voters' preferences

# a) Correlation between male and female preferences 
figure_1a <- ggscatter(data = parl_votes %>%
                      distinct(ref_id, .keep_all = TRUE), 
                    x = "vox_female_vote", y = "vox_male_vote", 
                    color = "vox_gender_majority_diff" ,shape = "vox_gender_majority_diff",
                    cor.coef = FALSE, xlab = "yes-share female electorate", ylab = "yes-share male electorate") + 
  geom_abline(intercept = 0, slope = 1, colour = "grey", linewidth = 12, alpha=0.3) +
  geom_abline(intercept = 0, slope = 1) +
  scale_color_manual(values = c("#96b6ff","#c882fa")) +
  scale_y_continuous(labels = percent, limits = c(0,1), expand = c(0, 0)) +
  scale_x_continuous(labels = percent, limits = c(0,1), expand = c(0, 0))+ 
  coord_fixed()+
  theme_minimal() +
  theme(legend.position = "none",
        text = element_text(size=14))
figure_1a

# b) Box plot differences in preferences

figure_1b <- ggplot(parl_votes %>%
                   distinct(ref_id, .keep_all = TRUE), 
                 aes(x = 0.5, y = vox_gender_diff_abs)) + 
  geom_boxplot() +
  geom_point(aes(x=yes_share_ch, y=vox_gender_diff_abs, color = vox_gender_majority_diff, shape = vox_gender_majority_diff))  +
  scale_color_manual(values = c("#96b6ff","#c882fa")) +
  stat_summary(fun=mean,geom="point", 
               shape=8, size=4, show.legend=FALSE) + 
  coord_flip() +
  scale_y_continuous(labels = function(x) paste0(x*100, "pp"), limits = c(0,0.3),expand = c(0, 0)) +
  scale_x_continuous(labels = percent, limits = c(0,1), expand = c(0, 0))+ 
  labs(y="absolute difference in yes-share between\nwomen and men in electorate",
       x="national yes-share") +
  theme_minimal() +
  theme(axis.ticks.y  = element_blank(),
        legend.position = "none",
        text = element_text(size=14))

figure_1b

#---------------------
# C) Empirical Analysis ####
#---------------------

#---------------------
# i) Table 1
#---------------------
# Effect of female voters' preferences on decisions by female and male legislators

# Specification (1)
i101 <- feglm(decision ~ 1
              + female*vox_female_vote
              ,
              data=parl_votes, family=binomial(link="logit"), vcov = ~mp_id)

# Specification (2)
i102 <- feglm(decision ~ 1
              + female*vox_female_vote
              + yes_share_canton 
              ,
              data=parl_votes, family=binomial(link="logit"), vcov = ~mp_id)

# Specification (3)
i103 <- feglm(decision ~ 1
              + female*vox_female_vote
              + party_yes 
              ,
              data=parl_votes, family=binomial(link="logit"), vcov = ~mp_id)

# Specification (4)
i104 <- feglm(decision ~ 1
              + female*vox_female_vote
              + yes_share_canton 
              + party_yes 
              ,
              data=parl_votes, family=binomial(link="logit"), vcov = ~mp_id)

# Specification (5)
i105 <- feglm(decision ~ 1
              + female*vox_female_vote
              + yes_share_canton 
              + party_yes 
              + age + age_squared + years_in_parliament + years_in_parliament_squared 
              + doctor_master
              + party_left + party_right
              ,
              data=parl_votes, family=binomial(link="logit"), vcov = ~mp_id)

# Specification (6)
i106 <- feglm(decision ~ 1
              + female*vox_female_vote
              + yes_share_canton
              + party_yes 
              + age + age_squared + years_in_parliament + years_in_parliament_squared 
              + doctor_master
              + party_left + party_right
              | canton_id + leg_period + ref_type
              ,
              data=parl_votes, family=binomial(link="logit"), vcov = ~mp_id)

i1est <- mget(ls(pattern = "i10"))

# Table 1
table_1 <- etable(i1est,
                 fitstat = ~ pr2 + n, digits = 3, digits.stats = 3, se.row = FALSE,
                 dict=c(yes_share_canton="Constituency Preferences",
                      vox_female_vote="Female Preferences",
                      party_yes = "Party Yes",
                      female ="Female",
                      canton_id="Canton",
                      leg_period="Period",
                      ref_type="Ref. type",
                      decision="MP YES"),
                 group = list(Controls = "age"),
                 keep = c("(Female)*(Female Preferences)","Constituency Preferences","Female Preferences","Female","Part Yes"))
table_1
#---------------------
# ii) Figure 2
#---------------------
# Discrete effect of an increase in female voters' preferences on female and male legislators' decisions

d1 <- i1est %>% 
  map(~as_tibble(discrete_triple(.))) %>% 
  set_names( nm = paste0("Model ", 1:length(i1est))) %>% 
  bind_rows(.id = "Model")
d1

figure_2 <- ggplot(data = d1 %>%
                 filter(effect != "difference"),
               aes(x = Model, y = Estimate, shape = effect, color = effect)) +
  geom_point(position = position_dodge(width = 0.5), size = 3) +
  geom_errorbar(aes(ymin = Estimate - SE * conf_multiplier, ymax = Estimate + SE * conf_multiplier),
                position = position_dodge(width = 0.5), width = 0.2, linewidth = 0.5) +
  theme_minimal() +
  labs(y = "discrete effect *Female Preferences*") +
  scale_y_continuous(label = percent) +
  scale_shape_manual(labels = c("female legislator","male legislator" ), values = c(17,19)) +
  scale_color_manual(labels = c("female legislator","male legislator" ),values = c("#c882fa", "#96b6ff")) +
  theme(legend.position = "bottom", 
        legend.title = element_blank(),
        axis.title.y = ggtext::element_markdown(),
        text=element_text(size=11),
        legend.text = element_text(size=11)) +
  xlab(NULL) 

figure_2

#---------------------
# iii) Table 2
#---------------------
# Effect of female voters' preferences on decisions by female and male legislators - social and gender-specific policies

# Create necessary subsets
parl_votes_social <- subset(parl_votes, top_10_social==1)
parl_votes_nonsocial <- subset(parl_votes, top_10_social==0)
parl_votes_divergent <- subset(parl_votes,  vox_gender_majority_diff==1)
parl_votes_nondivergent <- subset(parl_votes,  vox_gender_majority_diff==0)

# Specification (1)
i201 <- feglm(decision ~ 1
              + female*vox_female_vote
              + yes_share_canton
              + age +age_squared + years_in_parliament + years_in_parliament_squared
              + party_yes + doctor_master
              + party_left + party_right
              | canton_id + ref_type + leg_period
              ,
              data=parl_votes_social,family=binomial(link="logit"), vcov = ~mp_id)

# Specification (2)
i202 <- feglm(decision ~ 1
              + female*vox_female_vote
              + yes_share_canton
              + age +age_squared + years_in_parliament + years_in_parliament_squared
              + party_yes + doctor_master
              + party_left + party_right
              | canton_id + ref_type + leg_period
              ,
              data=parl_votes_nonsocial,family=binomial(link="logit"), vcov = ~mp_id)

# Specification (3)
i203 <- feglm(decision ~ 1
              + female*vox_female_vote*top_10_social
              + yes_share_canton
              + party_yes
              + age +age_squared + years_in_parliament + years_in_parliament_squared
              + doctor_master
              + party_left + party_right
              | canton_id + ref_type + leg_period
              ,
              data=parl_votes,family=binomial(link="logit"), vcov = ~mp_id)

# Specification (4)
i204 <- feglm(decision ~ 1
              + female*vox_female_vote
              + yes_share_canton
              + age +age_squared + years_in_parliament + years_in_parliament_squared
              + party_yes + doctor_master
              + party_left + party_right
              | canton_id + ref_type + leg_period
              ,
              data=parl_votes_divergent,family=binomial(link="logit"), vcov = ~mp_id)

# Specification (5)
i205 <- feglm(decision ~ 1
              + female*vox_female_vote
              + yes_share_canton
              + age +age_squared + years_in_parliament + years_in_parliament_squared
              + party_yes + doctor_master
              + party_left + party_right
              | canton_id + ref_type + leg_period
              ,
              data=parl_votes_nondivergent,family=binomial(link="logit"), vcov = ~mp_id)

i2est <- mget(ls(pattern = "i20"))

# Table 2
table_2 <- etable(i2est, 
       fitstat = ~ pr2 + n, digits = 3, digits.stats = 3, se.row = FALSE,
       dict=c(yes_share_canton="Constituency Preferences",
              vox_female_vote="Female Preferences", 
              female ="Female",
              vox_gender_majority_diff="Divided Preferences",
              party_yes = "Party Yes",
              canton_id="Canton",
              leg_period="Period",
              ref_type="Ref. type",
              decision="MP YES",
              top_10_social="Social topic"),
       group = list(Controls = "age"),
       headers = list(":_:"=list("Panel A: Social vs. non-social policy"=3, "Panel B: Divided gender preferences"=2)),
       keep = c("Constituency Preferences","Social topic","(Female)*(Female Preferences)*(Divided Preferences)",
                "(Female)*(Female Preferences)*(Socia lIssue)","Female Preferences","Female","Party Yes"))

table_2

# Discrete effects Female Preferences
d2 <- i2est %>% 
  map(~as_tibble(discrete_triple(.))) %>% 
  set_names( nm = paste0("Model ", 1:length(i2est))) %>% 
  bind_rows(.id = "Model")
d2


#---------------------
# iv) Figure 3
#---------------------
# Heterogeneity for social policy issues

# 1. Age
# Subset below Median
h101 <- feglm(decision ~ 1
              + female*vox_female_vote
              + yes_share_canton
              + age +age_squared + years_in_parliament + years_in_parliament_squared
              + doctor_master + party_yes
              + party_left + party_right
              |canton_id + leg_period + ref_type
              ,
              data=subset(parl_votes_social, age < median(parl_votes$age)),binomial(link="logit"), vcov = ~mp_id)

# Subset above Median
h102 <- feglm(decision ~ 1
              + female*vox_female_vote
              + yes_share_canton
              + age +age_squared + years_in_parliament + years_in_parliament_squared
              + doctor_master + party_yes
              + party_left + party_right
              |canton_id + leg_period + ref_type
              ,
              data=subset(parl_votes_social, age > median(parl_votes$age)), binomial(link="logit"), vcov = ~mp_id)

# 2. Vote on female suffrage in 1971 
# Suffrage accepted
h201 <- feglm(decision ~ 1
              + female*vox_female_vote
              + yes_share_canton
              + age +age_squared + years_in_parliament + years_in_parliament_squared
              + doctor_master + party_yes
              + party_left + party_right
              |canton_id + leg_period + ref_type
              ,
              data=subset(parl_votes_social,female_suffrage_accepted==1), binomial(link="logit"), vcov = ~mp_id)

# Suffrage rejected
h202 <- feglm(decision ~ 1
              + female*vox_female_vote
              + yes_share_canton
              + age +age_squared + years_in_parliament + years_in_parliament_squared
              + doctor_master + party_yes
              + party_left + party_right
              |canton_id + leg_period + ref_type
              ,
              data=subset(parl_votes_social, female_suffrage_accepted == 0), binomial(link="logit"), vcov = ~mp_id)


# 3. Share of women per party
# Subset Below median
h301 <- feglm(decision ~ 1
              + female*vox_female_vote
              + yes_share_canton
              + age +age_squared + years_in_parliament + years_in_parliament_squared
              + doctor_master + party_yes
              + party_left + party_right
              |canton_id + leg_period + ref_type
              ,
              data=subset(parl_votes_social, female_share_party <= mean(parl_votes$female_share_party)), binomial(link="logit"), vcov = ~mp_id)

# Subset Below median
h302 <- feglm(decision ~ 1
              + female*vox_female_vote
              + yes_share_canton
              + age +age_squared + years_in_parliament + years_in_parliament_squared
              + doctor_master + party_yes
              + party_left + party_right
              |canton_id + leg_period + ref_type
              ,
              data=subset(parl_votes_social, female_share_party > mean(parl_votes$female_share_party)), binomial(link="logit"), vcov = ~mp_id)


# 4. Electoral System
# Subset Proportional
h401 <- feglm(decision ~ 1
              + female*vox_female_vote
              + yes_share_canton
              + age +age_squared + years_in_parliament + years_in_parliament_squared
              + doctor_master + party_yes
              + party_left + party_right
              |canton_id + leg_period + ref_type
              ,
              data=subset(parl_votes_social, council_id == 0 | seats_canton >= 2), binomial(link="logit"), vcov = ~mp_id)

# Subset majoritarian
h402 <- feglm(decision ~ 1
              + female*vox_female_vote
              + yes_share_canton
              + age +age_squared + years_in_parliament + years_in_parliament_squared
              + doctor_master + party_yes
              + party_left + party_right
              |canton_id + leg_period + ref_type
              ,
              data=subset(parl_votes_social, council_id == 1 | seats_canton < 2), binomial(link="logit"), vcov = ~mp_id)



# Subplots
figure_3_a <- plot_summs(h101, h102, 
                     robust = TRUE,
                     ci_level = 0.95,
                     coefs = c("(Female)*(Female Preferences)"="female:vox_female_vote"),
                     legend.title = "Age in years",
                     point.shape = c(19,17),
                     point.size = 4,
                     colors = c("#96b6ff","#c882fa"),
                     line.size = c(0.5,1),
                     model.names = c("below median", "above median"))+
  scale_x_continuous(name = "Coefficient for *(Female)***(Female Preferences)*",limits=c(-1.5,4), n.breaks=6) +
  theme_minimal()+
  theme(axis.text.x = element_blank(),
        axis.title.x = element_blank(),
        legend.direction = "vertical",
        axis.title.y = ggtext::element_markdown(size = 12),
        legend.text = element_text(size = 11),
        panel.grid.major.x = element_blank())+
  coord_flip()

figure_3_b <- plot_summs(h201,h202,
                     robust = TRUE,
                     ci_level = 0.95,
                     coefs = c("(Female)*(Female Preferences)"="female:vox_female_vote"),
                     legend.title = "1971 vote on female suffrage",
                     point.shape = c(19,17),
                     point.size = 4,
                     line.size = c(0.5,1),
                     model.names = c("rejected", "accepted"),
                     colors = c("#96b6ff","#c882fa"))+
  scale_x_continuous(limits=c(-1.5,4), n.breaks=6) +
  theme_minimal()+
  theme(axis.text.y = element_blank(),
        axis.text.x = element_blank(),
        axis.title.x = element_blank(),
        axis.title.y = element_blank(),
        legend.text = element_text(size = 11),
        legend.direction = "vertical",
        panel.grid.major.x = element_blank())+
  coord_flip()

figure_3_c <- plot_summs(h301, h302,
                     robust = TRUE,
                     ci_level = 0.95,
                     coefs = c("(Female)*(Female Preferences)"="female:vox_female_vote"),
                     legend.title = "Female share in party",
                     point.shape = c(19,17),
                     point.size = 4,
                     line.size = c(0.5,1),
                     model.names = c("below mean", "above mean"),
                     colors = c("#96b6ff","#c882fa"))+
  scale_x_continuous(limits=c(-1.5,4), n.breaks=6) +
  theme_minimal()+
  theme(axis.text.y = element_blank(),
        axis.text.x = element_blank(),
        axis.title.y = element_blank(),
        legend.text = element_text(size = 11),
        axis.title.x = element_blank(),
        legend.direction = "vertical",
        panel.grid.major.x = element_blank())+
  coord_flip()

figure_3_d <- plot_summs(h401, h402,
                     robust = TRUE,
                     ci_level = 0.95,
                     coefs = c("(Female)*(Female Preferences)"="female:vox_female_vote"),
                     legend.title = "Electoral system",
                     point.shape = c(19,17),
                     point.size = 4,
                     line.size = c(0.5,1),
                     model.names = c("proportional", "majoritarian"),
                     colors = c("#96b6ff","#c882fa"))+
  scale_x_continuous(limits=c(-1.5,4), n.breaks=6) + 
  theme_minimal()+
  theme(axis.text.y = element_blank(),
        axis.text.x = element_blank(),
        legend.text = element_text(size = 11),
        axis.title.y = element_blank(),
        axis.title.x = element_blank(),
        legend.direction = "vertical",
        panel.grid.major.x = element_blank())+
  coord_flip()

# Figure 3
figure_3 <- ggarrange(figure_3_a, figure_3_b, figure_3_c, figure_3_d,
                  legend = "bottom",
                  ncol = 4, nrow = 1,
                  widths = c(1.2,1,1,1))
figure_3




###############################################################################
# Online Appendix
###############################################################################

#---------------------
# D) Additional Tables
#---------------------

#---------------------
# i) Table A1
#---------------------
# Descriptive statistics

table_a1 <- sumtable(parl_votes,
                   vars = c('decision',
                            'vox_female_vote',
                            'vox_male_vote',
                            'female_match',
                            'yes_share_canton',
                            'party_yes',
                            'female',
                            'age',
                            'doctor_master',
                            'years_in_parliament',
                            'council_id',
                            'party_left',
                            'party_center',
                            'party_right',
                            'female_suffrage_accepted',
                            'female_share_party',
                            'top_10_social'),
                   summ = c('notNA(x)', 'mean(x)', 'sd(x)'),
                   summ.names =c ('N', 'Mean', 'Sd'),
                   out = 'return',
                   digits = 3,
                   factor.numeric = TRUE)
table_a1

#---------------------
# ii) Table A2
#---------------------
# Female and male voters' preferences for each referendum

table_a2 <- parl_votes %>%
  distinct(ref_id, .keep_all = TRUE) %>%
  select(ref_id, vox_female_vote, vox_male_vote) %>%
  mutate(vox_female_vote = vox_female_vote*100,
         vox_male_vote = vox_male_vote*100,
         diff = vox_female_vote-vox_male_vote)
print(table_a2, n = Inf)

#---------------------
# iii) Table A3
#---------------------
# Effect of female voters' preferences on decisions by female and male legislators - 
# Controlling for the electorate's preferences at the national level, as measured by referendum results

# Specification (1)
a301 <- feglm(decision ~ 1
              + female*vox_female_vote
              ,
              data=parl_votes, family=binomial(link="logit"), vcov = ~mp_id)

# Specification (2)
a302 <- feglm(decision ~ 1
              + female*vox_female_vote
              + yes_share_ch 
              ,
              data=parl_votes, family=binomial(link="logit"), vcov = ~mp_id)
# Specification (3)
a303 <- feglm(decision ~ 1
              + female*vox_female_vote
              + party_yes 
              ,
              data=parl_votes, family=binomial(link="logit"), vcov = ~mp_id)

# Specification (4)
a304 <- feglm(decision ~ 1
              + female*vox_female_vote
              + yes_share_ch 
              + party_yes 
              ,
              data=parl_votes, family=binomial(link="logit"), vcov = ~mp_id)

# Specification (5)
a305 <- feglm(decision ~ 1
              + female*vox_female_vote
              + yes_share_ch 
              + party_yes 
              + age + age_squared + years_in_parliament + years_in_parliament_squared 
              + doctor_master
              + party_left + party_right
              ,
              data=parl_votes, family=binomial(link="logit"), vcov = ~mp_id)

# Specification (6)
a306 <- feglm(decision ~ 1
              + female*vox_female_vote
              + yes_share_ch
              + party_yes 
              + age + age_squared + years_in_parliament + years_in_parliament_squared 
              + doctor_master
              + party_left + party_right
              | canton_id + leg_period + ref_type
              ,
              data=parl_votes, family=binomial(link="logit"), vcov = ~mp_id)

# Table A3
table_a3 <- etable(a301, a302, a303, a304, a305,a306,
                 fitstat = ~ pr2 + n, 
                 digits = 3,
                 digits.stats = 3,
                 se.row = FALSE,
                 dict=c(yes_share_ch="CH Preferences",
                        vox_female_vote="Female Preferences",
                        party_yes = "Party Yes",
                        female ="Female",
                        canton_id="Canton",
                        leg_period="Period",
                        ref_type="Ref. type",
                        decision="MPYES"),
                 group = list(Controls = "age"),
                 keep = c("(Female)*(Female Preferences)","CH Preferences","Female Preferences","Female","Party Yes"))
table_a3


#---------------------
# iv) Table A4
#---------------------
# Effect of female preferences on decisions by female and male legislators - 
# Controlling for the electorate's preferences at the national level, as measured by postreferendum surveys

# Specification (1)
a401 <- feglm(decision ~ 1
              + female*vox_female_vote
              ,
              data=parl_votes, family=binomial(link="logit"), vcov = ~mp_id)

# Specification (2)
a402 <- feglm(decision ~ 1
              + female*vox_female_vote
              + vox_yes_share 
              ,
              data=parl_votes, family=binomial(link="logit"), vcov = ~mp_id)
# Specification (3)
a403 <- feglm(decision ~ 1
              + female*vox_female_vote
              + party_yes 
              ,
              data=parl_votes, family=binomial(link="logit"), vcov = ~mp_id)

# Specification (4)
a404 <- feglm(decision ~ 1
              + female*vox_female_vote
              + vox_yes_share 
              + party_yes 
              ,
              data=parl_votes, family=binomial(link="logit"), vcov = ~mp_id)

# Specification (5)
a405 <- feglm(decision ~ 1
              + female*vox_female_vote
              + vox_yes_share 
              + party_yes 
              + age + age_squared + years_in_parliament + years_in_parliament_squared 
              + doctor_master
              + party_left + party_right
              ,
              data=parl_votes, family=binomial(link="logit"), vcov = ~mp_id)

# Specification (6)
a406 <- feglm(decision ~ 1
              + female*vox_female_vote
              + vox_yes_share
              + party_yes 
              + age + age_squared + years_in_parliament + years_in_parliament_squared 
              + doctor_master
              + party_left + party_right
              | canton_id + leg_period + ref_type
              ,
              data=parl_votes, family=binomial(link="logit"), vcov = ~mp_id)

# Table A4
table_a4 <- etable(a401, a402, a403, a404, a405,a406,
                 fitstat = ~ pr2 + n, 
                 digits = 3,
                 digits.stats = 3,
                 se.row = FALSE,
                 dict=c(vox_yes_share="CH Preferences VOX",
                        vox_female_vote="Female Preferences",
                        party_yes = "Party Yes",
                        female ="Female",
                        canton_id="Canton",
                        leg_period="Period",
                        ref_type="Ref. type",
                        decision="MP YES"),
                 group = list(Controls = "age"),
                 keep = c("(Female)*(Female Preferences)","CH Preferences VOX","Female Preferences","Female","Party Yes"))
table_a4

#---------------------
# v) Table A5
#---------------------
# Heterogeneous effect of female voters' preferences on female and male legislators' social policy decisions

table_a5 <- etable(h101, h102, h201, h202, h301, h302, h401, h402, 
                 fitstat = ~ pr2 + n, 
                 digits = 3,
                 digits.stats = 3,
                 se.row = FALSE,
                 dict=c(yes_share_canton="Constituency Preferences",
                        vox_female_vote="Female Preferences", 
                        female ="Female",
                        canton_id="Canton",
                        party_yes = "Party Yes",
                        leg_period="Period",
                        ref_type="Ref. type",
                        decision="MPYES"),
                 group = list(Controls = "years_in_parliament"),
                 keep = c("Constituency Preferences","Female Preferences","Female","Party Yes"))
table_a5

#---------------------
# vi) Table A6
#---------------------
# Heterogeneous effect of female voters' preferences on female and male legislators' decisions on non-social policy proposals

# 1. Age
# Subset below Median
a601 <- feglm(decision ~ 1
              + female*vox_female_vote
              + yes_share_canton
              + age +age_squared + years_in_parliament + years_in_parliament_squared
              + doctor_master + party_yes
              + party_left + party_right
              |canton_id + leg_period + ref_type
              ,
              data=subset(parl_votes_nonsocial, age < median(parl_votes$age)),binomial(link="logit"), vcov = ~mp_id)


# Subset above Median
a602 <- feglm(decision ~ 1
              + female*vox_female_vote
              + yes_share_canton
              + age +age_squared + years_in_parliament + years_in_parliament_squared
              + doctor_master + party_yes
              + party_left + party_right
              |canton_id + leg_period + ref_type
              ,
              data=subset(parl_votes_nonsocial, age > median(parl_votes$age)), binomial(link="logit"), vcov = ~mp_id)

# 2. Vote on female suffrage in 1971 
# Suffrage accepted
a603 <- feglm(decision ~ 1
              + female*vox_female_vote
              + yes_share_canton
              + age +age_squared + years_in_parliament + years_in_parliament_squared
              + doctor_master + party_yes
              + party_left + party_right
              |canton_id + leg_period + ref_type
              ,
              data=subset(parl_votes_nonsocial,female_suffrage_accepted==1),binomial(link="logit"), vcov = ~mp_id)

# Sufferage rejected
a604 <- feglm(decision ~ 1
              + female*vox_female_vote
              + yes_share_canton
              + age +age_squared + years_in_parliament + years_in_parliament_squared
              + doctor_master + party_yes
              + party_left + party_right
              |canton_id + leg_period + ref_type
              ,
              data=subset(parl_votes_nonsocial,female_suffrage_accepted==0),binomial(link="logit"), vcov = ~mp_id)

# 3. Share of women per party
# Subset Below median
a605 <- feglm(decision ~ 1
              + female*vox_female_vote
              + yes_share_canton
              + age +age_squared + years_in_parliament + years_in_parliament_squared
              + doctor_master + party_yes
              + party_left + party_right
              |canton_id + leg_period + ref_type
              ,
              data=subset(parl_votes_nonsocial, female_share_party <= mean(parl_votes$female_share_party)),binomial(link="logit"), vcov = ~mp_id)

# Subset Below median
a606 <- feglm(decision ~ 1
              + female*vox_female_vote
              + yes_share_canton
              + age +age_squared + years_in_parliament + years_in_parliament_squared
              + doctor_master + party_yes
              + party_left + party_right
              |canton_id + leg_period + ref_type
              ,
              data=subset(parl_votes_nonsocial, female_share_party > mean(parl_votes$female_share_party)),binomial(link="logit"), vcov = ~mp_id)

# 4. Electoral System
# Subset Proportional
a607 <- feglm(decision ~ 1
              + female*vox_female_vote
              + yes_share_canton
              + age +age_squared + years_in_parliament + years_in_parliament_squared
              + doctor_master + party_yes
              + party_left + party_right
              |canton_id + leg_period + ref_type
              ,
              data=subset(parl_votes_nonsocial,council_id==0 | seats_canton >=2),binomial(link="logit"), vcov = ~mp_id)


# Subset majoritarian
a608 <- feglm(decision ~ 1
              + female*vox_female_vote
              + yes_share_canton
              + age +age_squared + years_in_parliament + years_in_parliament_squared
              + doctor_master + party_yes
              + party_left + party_right
              |canton_id + leg_period + ref_type
              ,
              data=subset(parl_votes_nonsocial,council_id ==1| seats_canton <2),binomial(link="logit"), vcov = ~mp_id)

# Table A6
table_a6 <- etable(a601, a602, a603, a604, a605, a606, a607, a608, 
                 fitstat = ~ pr2 + n, 
                 digits = 3,
                 digits.stats = 3,
                 se.row = FALSE,
                 dict=c(yes_share_canton="Constituency Preferences",
                        vox_female_vote="Female Preferences", 
                        female ="Female",
                        canton_id="Canton",
                        party_yes = "Party Yes",
                        leg_period="Period",
                        ref_type="Ref. type",
                        decision="MPYES"),
                 group = list(Controls = "years_in_parliament"),
                 keep = c("Constituency Preferences","Female Preferences","Female","Party Yes"))
table_a6

#---------------------
# vii) Table A7
#---------------------
# Heterogeneous effect of female voters' preferences on female and male legislators' social policy decisions - Triple interaction

a701 <- feglm(decision ~ 1
              + female*vox_female_vote*age_high
              + yes_share_canton
              + age +age_squared + years_in_parliament + years_in_parliament_squared
              + doctor_master + party_yes
              + party_left + party_right
              |canton_id + leg_period + ref_type
              ,
              data=subset(parl_votes_social),binomial(link="logit"))
a702 <- feglm(decision ~ 1
              + female*vox_female_vote*age
              + yes_share_canton
              + age_squared + years_in_parliament + years_in_parliament_squared
              + doctor_master + party_yes
              + party_left + party_right
              |canton_id + leg_period + ref_type
              ,
              data=subset(parl_votes_social),binomial(link="logit"))
a703 <- feglm(decision ~ 1
              + female*vox_female_vote*female_suffrage_accepted
              + yes_share_canton
              + age +age_squared + years_in_parliament + years_in_parliament_squared
              + doctor_master + party_yes
              + party_left + party_right
              |canton_id + leg_period + ref_type
              ,
              data=subset(parl_votes_social),binomial(link="logit"))
a704 <- feglm(decision ~ 1
              + female*vox_female_vote*proportional
              + yes_share_canton
              + age +age_squared + years_in_parliament + years_in_parliament_squared
              + doctor_master + party_yes
              + party_left + party_right
              |canton_id + leg_period + ref_type
              ,
              data=subset(parl_votes_social),binomial(link="logit"))
a705 <- feglm(decision ~ 1
              + female*vox_female_vote*female_share_party_high
              + yes_share_canton
              + age +age_squared + years_in_parliament + years_in_parliament_squared
              + doctor_master + party_yes
              + party_left + party_right
              |canton_id + leg_period + ref_type
              ,
              data=subset(parl_votes_social),binomial(link="logit"))
a706 <- feglm(decision ~ 1
              + female*vox_female_vote*female_share_party
              + yes_share_canton
              + age +age_squared + years_in_parliament + years_in_parliament_squared
              + doctor_master + party_yes
              + party_left + party_right
              |canton_id + leg_period + ref_type
              ,
              data=subset(parl_votes_social),binomial(link="logit"))


table_a7 <- etable(a701, a702, a703,  a704, a705, a706,
                 fitstat = ~ pr2 + n, 
                 digits = 3,
                 digits.stats = 3,
                 se.row = FALSE,
                 dict=c(yes_share_canton="Constituency Preferences",
                        vox_female_vote="Female Preferences", 
                        female ="Female",
                        canton_id="Canton",
                        party_yes = "Party Yes",
                        leg_period="Period",
                        ref_type="Ref. type",
                        decision="MPYES"),
                 group = list(Controls = "years_in_parliament"),
                 keep = c("Constituency Preferences","Female Preferences","Female","Party Yes"))
table_a7

#---------------------
# viii) Table A8
#---------------------
# Effect of gender on congruence with the majority of female voters

# Specification (1)
a801 <- feglm(female_match ~ 1
              + female
              + age + age_squared + years_in_parliament + years_in_parliament_squared + doctor_master
              |canton_id + leg_period + ref_type
              ,
              data=parl_votes, family=binomial(link="logit"))

# Specification (2)
a802 <- feglm(female_match ~ 1
              + female
              + age + age_squared + years_in_parliament + years_in_parliament_squared + doctor_master
              + party_left + party_right
              |canton_id + leg_period + ref_type
              ,
              data=parl_votes, family=binomial(link="logit"))

# Specification (3)
a803 <- feglm(female_match ~ 1
              + female
              + age + age_squared + years_in_parliament + years_in_parliament_squared + doctor_master
              |canton_id + leg_period + ref_type
              ,
              data=parl_votes_social, family=binomial(link="logit")) 

# Specification (4)
a804 <- feglm(female_match ~ 1
              + female
              + age + age_squared + years_in_parliament + years_in_parliament_squared + doctor_master
              + party_left + party_right
              |canton_id + leg_period + ref_type
              ,
              data=parl_votes_social, family=binomial(link="logit"))

# Specification (5)
a805 <- feglm(female_match ~ 1
              + female
              + age + age_squared + years_in_parliament + years_in_parliament_squared + doctor_master
              |canton_id + leg_period + ref_type
              ,
              data=subset(parl_votes, vox_gender_majority_diff ==1 ), family=binomial(link="logit"))

# Specification (6)
a806 <- feglm(female_match ~ 1
              + female
              + age + age_squared + years_in_parliament + years_in_parliament_squared + doctor_master
              + party_left + party_right
              |canton_id + leg_period + ref_type
              ,
              data=subset(parl_votes, vox_gender_majority_diff ==1 ), family=binomial(link="logit"))


table_a8 <- etable(a801,a802,a803,a804,a805,a806,
                 fitstat = ~ pr2 + n, 
                 digits = 3,
                 digits.stats = 3,
                 se.row = FALSE,
                 dict=c(canton_match="MP = Constituency",
                        female_match = "MP = Female",
                        male_match = "MP = Male",
                        female ="Female",
                        canton_id="Canton",
                        party_yes = "Party Yes",
                        party_left = "PartyLeft",
                        party_right = "PartyRight",
                        leg_period="Period",
                        ref_type="Ref. type",
                        decision="MPYES"),
                 group = list(MP_Controls = "years_in_parliament"),
                 keep = c("Female","PartyLeft","PartyRight","Constituency = Country","Female electorate = Country", "Male electorate = Country"))
table_a8

#---------------------
# E) Additional Figures
#---------------------

#---------------------
# i) Figure B1
#---------------------
# Correlation between female and constituency preferences

# Create dataset dith distinct referendum-canton pairs
unique_ref_constituency_id <- parl_votes %>% 
  distinct(ref_id,canton_id, .keep_all = TRUE) %>%
  mutate(f_cant_diff_abs = abs(vox_female_vote-yes_share_canton))

figure_b1 <- ggscatter(data=unique_ref_constituency_id, x = "vox_female_vote", y = "yes_share_canton", color = "f_cant_diff_abs",size = 1,
                    cor.coef = FALSE, 
                    xlab = "yes-share female electorate", ylab = "yes-share constituency") + 
  geom_abline(intercept = 0, slope = 1, colour = "grey", linewidth = 14, alpha=0.3) +
  geom_abline(intercept = 0, slope = 1) +
  scale_y_continuous(labels = percent, limits = c(0,1),expand = c(0, 0)) +
  scale_x_continuous(labels = percent, limits = c(0,1),expand = c(0, 0))+ 
  scale_color_gradient(low = "#00A087B2", high = "red",guide = "none") +
  coord_fixed()+
  theme_minimal() +
  theme(text = element_text(size=12))
figure_b1


#---------------------
# ii) Figure B2
#---------------------
# Heterogeneity and robustness - Non-social policies

# Subplots
# Age
figure_b2_a <- plot_summs(a601, a602,
                       robust = TRUE,
                       ci_level = 0.95,
                       coefs = c("(Female)*(Female Preferences)"="female:vox_female_vote"),
                       legend.title = "Age in years",
                       line.size = c(0.5,1),
                       point.shape = c(19,17),
                       point.size = 5,
                       colors = c("#96b6ff","#c882fa"),
                       model.names = c("below median", "above median"))+
  scale_x_continuous(name = "Coefficient for *(Female)***(Female Preferences)*",limits=c(-3.5,3.5), n.breaks=8) +
  theme_minimal()+
  theme(axis.text.x = element_blank(),
        axis.title.x = element_blank(),
        legend.direction = "vertical",
        axis.title.y = ggtext::element_markdown(size = 12),
        legend.text = element_text(size = 11),
        legend.position = "top",
        panel.grid.major.x = element_blank())+
  coord_flip()

# Female suffrage
figure_b2_b <- plot_summs(a603,a604,
                       robust = TRUE,
                       ci_level = 0.95,
                       coefs = c("(Female)*(Female Preferences)"="female:vox_female_vote"),
                       legend.title = "1971 vote on female suffrage",
                       line.size = c(0.5,1),
                       point.shape = c(19,17),
                       point.size = 5,
                       colors = c("#96b6ff","#c882fa"),
                       model.names = c("accepted", "rejected")) +
  scale_x_continuous(limits=c(-3.5,3.5), n.breaks=8) +
  theme_minimal()+
  theme(axis.text.y = element_blank(),
        axis.text.x = element_blank(),
        axis.title.x = element_blank(),
        axis.title.y = element_blank(),
        legend.text = element_text(size = 11),
        legend.direction = "vertical",
        legend.position = "top",
        panel.grid.major.x = element_blank())+
  coord_flip()

# Female share in party
figure_b2_c <- plot_summs(a605, a606,
                       robust = TRUE,
                       ci_level = 0.95,
                       coefs = c("(Female)*(Female Preferences)"="female:vox_female_vote"),
                       legend.title = "Female share in party",
                       line.size = c(0.5,1),
                       point.shape = c(19,17),
                       point.size = 5,
                       colors = c("#96b6ff","#c882fa"),
                       model.names = c("below mean", "above mean")) +
  scale_x_continuous(limits=c(-3.5,3.5), n.breaks=8) +
  theme_minimal()+
  theme(axis.text.y = element_blank(),
        axis.text.x = element_blank(),
        axis.title.y = element_blank(),
        legend.text = element_text(size = 11),
        axis.title.x = element_blank(),
        legend.direction = "vertical",
        legend.position = "top",
        panel.grid.major.x = element_blank())+
  coord_flip()

# Electoral rule
figure_b2_d <- plot_summs(a607, a608,
                       robust = TRUE,
                       ci_level = 0.95,
                       coefs = c("(Female)*(Female Preferences)"="female:vox_female_vote"),
                       legend.title = "Electoral system",
                       line.size = c(0.5,1),
                       point.shape = c(19,17),
                       point.size = 5,
                       colors = c("#96b6ff","#c882fa"),
                       model.names = c("proportional", "majoritarian")) +
  scale_x_continuous(limits=c(-3.5,3.5), n.breaks=8) +
  theme_minimal()+
  theme(axis.text.y = element_blank(),
        axis.text.x = element_blank(),
        legend.text = element_text(size = 11),
        axis.title.y = element_blank(),
        axis.title.x = element_blank(),
        legend.direction = "vertical",
        legend.position = "top",
        panel.grid.major.x = element_blank())+
  coord_flip()

# Figure B2
figure_b2 <- ggarrange(figure_b2_a, figure_b2_b, figure_b2_c, figure_b2_d,
                    legend = "bottom",
                    ncol = 4, nrow = 1,
                    widths = c(1.2,1,1,1))
figure_b2



